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ABSTRACT 

We investigate the effect of heating by luminosity sources in a simulation of clustered star formation. 
Our heating method involves a simplified continuum radiative transfer method that calculates the dust 
temperature. The gas temperature is set by the dust temperature. We present the results of four 
simulations, two simulations assume an isothermal equation of state and the two other simulations 
include dust heating. We investigate two mass regimes, i.e., 84M Q and 671M Q , using these two 
different energetics algorithms. The mass functions for the isothermal simulations and simulations 
which include dust heating are drastically different. In the isothermal simulation, we do not form 
any objects with masses above IMq. However, the simulation with dust heating, while missing some 
of the low-mass objects, forms high-mass objects (~ 20M Q ) which have a distribution similar to the 
Salpeter IMF. The envelope density profiles around the stars formed in our simulation match observed 
values around isolated, low-mass star-forming cores. We find the accretion rates to be highly variable 
and, on average, increasing with final stellar mass. By including radiative feedback from stars in a 
cluster-scale simulation, we have determined that it is a very important effect which drastically affects 
the mass function and yields important insights into the formation of massive stars. 
Subject headings: hydrodynamics — ISM: clouds — ISM: dust — methods: numerical — stars: for- 
mation 



1. INTRODUCTION 

One of the great puzzles in the field of star formation is 
the universality of the initi al mass function (IMF) in the 
Gala c tic environment (s ee Safoctcr 1955: Mill er fc Scald 
[19791: [Scaloin98^ : iKrounal j2t)T)51Chabrieri l2003h . specif- 
ically the slope at high mass and the characteristic 
turnover mass. One of the proposed explanations for the 
IMF's universality is similar thermal properties of the 
gas in the Galaxy in star-forming regions. The equation 
of state of the gas set by the thermal physics determines 
the level of fragmentat ion in a cloud and thus the mass 
function (|Larso 

The initial state of the gas that is likely to form stars 
is typically assumed to be dense, molecular gas which is 
isothermal, maintaining a temperature of ~10K. Simula- 
tions with isothermal equations of s tate have been used 
in an attempt to r ecreate the IMF (|Klessen etaD 11993 
iMartel et all l2006t hereafter MES06). These attempts 
have met with limited success due to the behavior of 
the Jeans mass, the basic unit of star formation, under 
the conditions of an isothermal equation of state. The 
Jeans mass is proportional to T 3 / 2 n -1 / 2 , where T is the 
gas temperature (K) and n is the gas number density 
(cm -3 ). For a constant temperature, as in the isothermal 
case, the Jeans mass decreases as the density increases, 
leading to perpetual fragmentation for a collapsing gas 
cloud. In reality, the fragmentation is believed to stop 
when the gas becomes optically thick as the density in- 
creases. This prevents radiation from escaping and cool- 
ing the gas back to the isothermal temperature. 

1 Department of Astronomy, University of Texas, Austin, TX 
78712 

2 Jet Propulsion Laboratory, California Institute of Technology, 
Pasadena, CA 91109 

3 Departement de physique, genie physique et optique, Univer- 
sity Laval, Quebec, QC, G1K 7P4, Canada 

4 Centre de Recherche en Astrophysique du Quebec 



In order to model this effect in simulations, a poly- 
tropic equation of state has been assumed, i.e. F oc p 7 . 
At low densities, 7 < 1 and at higher densities, 7 > 1. 
The density at which this transition occurs an d the values 
of 7 h ave been studied by many groups (see | Bate et all 
20031 iLi et all I200a iB atc 200 5Y I Jappsen et al.l 1200a 
Larsod[2005llBonnell et all 120^(1 IClark et alJl20Ciah . 

However, all of the studies mentioned above ignore a 
key ingredient in the attempt to re-create the IMF. Their 
assumption that the equation of state depends solely 
on density is likely to be valid only when the first gas 
cores condense out of the cloud and begin to collapse. 
Once these cores begin to generate their own energy via 
collapse, accretion, and deuterium/hydrogen burning, a 
spatially-uniform equation of state is no longer valid. 
The newly formed stars will heat the gas in the cloud 
and influence subsequent, nearby star formation. 

In order to include the energy from the forming 
stars in simulat i ons, r adiative transfer must be used. 
iKrumholz et~aP (|2007fl have run simulations in which 
they include the radiative energy from young stars. They 
include this effect in order to understand the formation 
of massive stars. Therefore, they study a small size scale 
(collapse from ~ 0.1 pc to ~ 10 AU scales) at very high 
densities (n ~ 10 9 cm~ 3 ) and form < 10 stars. They 
use the method of flux-limited diffusion (FLD) to cal- 
culate the effect of radiation on the surrounding dust. 
Then they assume the dust and gas are well-coupled 
to calculate the gas temperature. (This assumption is 
only valid in the high de nsity regions that they study 
(iGoldreich fc Kwanlll974f) ). As discussed in their paper, 
IKrumholz et al.l (|2007l ) assume gray radiation which will 
underestimate the true dust temperature. Also, their 
method of FLD is only accurate in very optically thick 
regimes. Despite these limitations, their method of radia- 
tive transfer is suitable for their study of the formation 
of massive stars. 
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iBatd (|2009bl) also includes FLD in his simulations. He 
forms less than 20 stars in his simulations which follow 
the collapse from scales of ~ O.lpc to 0.5 AU and the 
largest object formed has M < 2Mq. Because of the 
low mass of the objects formed in his simulations, he 
claims that it is valid to ignore the intrinsic stellar lu- 
minosity in his calculation of heating. He only considers 
energy generation via work on the gas. This only ac- 
counts for any accretion luminosity generated outside of 
the accretion radius. However, since most of the accre- 
tion luminosity is released when material accretes onto 
the star, his calculation of the accretion luminosity is an 
extreme underestimate of the total value. Hence, his re- 
sults only demonstrate the minimum possible effect of 
radiative feedback; 

Recent work by [Offner ct al. (2009) uses a method sim- 
ilar to the one used in lKrumholz et al.1 (120071). but studies 
a region similar in scale to that of lBatd ( 2009bD . They 
also form few objects (~ 15). They model a larger re- 
gion with a lower initial densit y. Their simulation differs 
from the work of IBatd (|2009b[ ) because they include nu- 
clear burni ng and accre tion luminosity (both effectively 
ignored bv IBatd (|2009bf )) which they state are the main 
sources of heating in their simulat ion. Therefore, they 
agree that the work of [Bate (20091]) only shows the min- 
imum possible effect of radiation on the star formation 
process. 

In our work, we attempt to model the effect of the heat- 
ing of dust and gas by young stars on the form of the IMF 
in a clustered environment. We study a cluster forming 
region with scales an order of magn i tude l arger than pre- 
viously studied bv lKrumholz et al.1 (|2007[ ) (~ 100 AU to 
1 pc). We are interested in the effect of our more realistic 
treatment of the temperature evolution on the mass func- 
tion. Are previous works which ignore the local effect of 
forming stars realistic? Or is their work only applicable 
to the very earliest stages of star formation when starless 
cores are forming? In order to answer these questions, 
we simulate the evolution of a star-forming region with 
a hydrodynamics + gravity code and allow the sources 
which form within our region to "turn-on" and heat the 
surrounding material. 
The form of r a diatiy e transfer that we use is described 
Urban et al.l (l2009f). It differs from the form used 
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Simulation Parameters 



Krumho lzet all (|2007t ) and [Bate] (|2009bh in that our 



in 
by 

method does not assume gray radiation, is applicable to 
a range of optical depths, and assumes that the matter 
distribution around each source is spherically symmetric. 
This latter approximation enables us to study a larger 
parameter space, which is necessary for modeling a clus- 
tered environment which is larger and where densities 
and optical depths are lower. The method we use cal- 
culates the dust temperature from the source luminosity 
and the density distribution around the source, using a 
grid of models gene rated by the spherical radiative trans- 
port co de DUSTY dNenkova et al.l l2000h. 

As in lKrumholz et all (|2007[ L we assume that the gas 
and dust are effectively coupled and that dust heating 
controls the temperature in the region. By ignoring other 
heating and cooling processes (such as cosmic rays and 
molecular cooling) , we can examine in detail the effect of 
adding local heating that can not be described using an 
analytic approximation. Our simulation with only dust 
heating provides a standard for future work. We expect 
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that including other cooling and heating mechanisms will 
only decrease the effect on the IMF that we will see in our 
simulation. Therefore the results we show in this paper 
represent the case of dust heating as the dominant term 
in the energy equation at all densities and temperatures. 

The simulations we present here are based on the work 
by MES06. We run a total of four new simulations. For 
all of our simulations, we assume that the gas is isother- 
mal until cores form at high densities. Then we consider 
two different energy transfer methods. In two of our 
simulations we assume that the gas is isothermal for the 
entire simulation runtime. In two different simulations, 
when the gas is dense and the first sink particle forms, 
we discard the assumption of isothermality and use the 
dust t emperature calculation described in lUrban et alJ 
(2009) to calculate the dust temperature and then set 
the gas temperature equal to the dust temperature. For 
each isothermal simulation and each simulation with dust 
heating, we have two size/mass scales, small and large, 
determined by the number of generations of particle split- 
ting. We refer to these small and large simulations by 
N gcn = 1 or 2, respectively 

In §21 we discuss the algorithms we use to solve the 
fluid equations and calculate the density profile, lumi- 
nosity, mass accretion rate, and dust temperature. In 
fj3j we describe the details of the simulations. In §4] we 
describe our results. This is followed by a discussion in 
§3 and our conclusions in §2) 

2. THE NUMERICAL ALGORITHM 

We use the smoothed particle hydrodynamics (SPH) 
algorithm with particle splitting and sink particles that 
is described in MES06. We have modified this algorithm 
to include the effect of the luminosity from forming stars 
on the dust/gas temperature in our simulation (see §2.31 
and gHH) . 

Other differences between this work and MES06 are 
the initial temperatures, 5K vs. 10K, and the threshold 
density contrasts for sink formation, 5,942 vs. 40,000, 
for this work and MES06, respectively. In the case of the 
isothermal simulations presented in MES06, the system 
modeled was scale-free. The parameters which deter- 
mined the properties of the simulation were the temper- 
ature and density contrast. Choosing an initial density to 
interpret the results then set the scale of the simulation. 
By this method, the authors were able to scale their sim- 
ulation to higher densities. In our work, we fix the initial 
conditions of temperature and density for all simulations. 
Therefore, the inclusion of particle splitting, used to in- 
crease the resolution of the simulations in MES06, will 
instead be used to study larger, more massive regions. 
This is because our simulation does not increase the den- 
sity resolution as in the case of the isothermal simulation 
in MES06, but instead increases the size of our simula- 
tion box while the density resolution is held fixed for 
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all simulations. For our small simulation {N gcn = 1), 
we simulate a region of volume ~ 0.5pc 3 with a mass 
of ~ 84A/0. For our large simulation, (iVgen = 2), the 
simulated region has a volume of ~ lpc 3 and a mass of 
~ 670M Q . Tabic [T] illustrates the difference in size and 
mass of our simulations with increasing levels of particle 
splitting. 

2.1. SPH 

We use a standard SPH algorithm (e.g., Monaghan 
I1992L and references therein), to simulate the evolution 
of a molecular cloud inside a cubic volume with periodic 
boundary conditions. This algorithm was modified to 
include particle splitting and sink particles. The Jeans 
criterion requires that a Jeans mass contains a certain 
minimum number of SPH particles to be properly re- 
solved, and prevent a sp urious numerical effect known as 
artificial fra gmentation dTrueloye et al.lll997t lBosslll998l : 
see however iHubber et al.l [2006). This translates into a 
condition between the mass, density, and specific inter- 
nal energy of particles. Whenever a particle violates this 
condition, the code splits the particle, replacing it by 
8 eq ual-mass particles located at the vertices of a cube 
(see iKitsionas fc Whitworthl 120021 MES06). Split parti- 
cles can later re-split if the condition is violated again. 
We allow for a maximum of N gen generations of split- 
ting. Hence, the mass ratio between the most and least 
massive particles is 8^™. 

Sink particles (or sinks) are created when the gas den- 
sity exceeds a density threshold p c . A group of particles, 
whose total mass equals one Jeans mass, are replaced 
by a single, massive sink particle, which has the abil- 
ity to grow by accreting gas particles. Boundary condi- 
tions are imposed at the interface between the sink and 
the surrounding gas. In MES06, we used the boundary 
conditions described in IBate et al.l ([l995). For this pa- 
pe r, we switched to th e boundary conditions described 
in IBromm et al.l ((2002). We should point out that the 
particular choice of boundary conditions is not critical 
in our simulations, because infall of gas onto sinks tends 
to be radia l and supersonic, making boundary conditions 
irrelevant (jBate et al.l fT995). It was easier to implement 
the dust physics (described in §2.4p i nto the algo r ithm if 
we use the boundary conditions of IBromm et al.l (2002), 
which is why we made th e switch. For details , we refer 
the reader to MES06 and IBromm et al.l (|2002f ). 

The sink particles in our simulation represent a star- 
forming core which may form a single star or a group of 
stars. Fragmentation within a sink particle due to pro- 
cesses that may occur at higher temperatures and den- 
sities are not modeled in our simulation. Therefore the 
masses of the sinks that we discuss throughout this paper 
only directly correspond to individual stars if no further 
fragmentation occurs within a sink. 

2.2. Density Profile 

In order to determine the dust temperature in our sim- 
ulation (as described later in §2.4ft , we must first calculate 
the density profile by a spherical average around individ- 
ual sinks. To do this, we bin the particles around each 
sink into concentric shells which hold exactly 25 gas par- 
ticles each, and calculate the density within each shell. 

The outer edge of the density distribution is set such 
that the mean density inside the outer radius is 200 times 



the initial background density. (This method is com- 
monly used in similar cosmological simulations, where 
this outer radius is referr ed to as the virial radius; see, 
e.g., iNavarro et al.lll995l .) We assume that it takes a 
minimum of 200 gas particles inside the outer radius to 
accurately determine the density profile. If this condition 
is not met, then we cannot calculate the density profile 
(we discuss the effect of this in £J2 .4|) . We show exam- 
ples of the density profile evolution in §4.21 The density 
profile is parameterized by n Q and a, 

n = n ° ( iooUau) %m ~ 3 - (1) 

Throughout this paper we will use the term, n, to repre- 
sent the number density of all particles (n = njj 2 + nn e ) 
assuming nn 2 /nn c = 5, which gives \l = 2.33 (this relates 
the number and volume density discussed in ij3.1[) . 

2.3. Luminosity and Mass Accretion Rate 

In order to calculate the luminosity of a sink par- 
ticle, we assume that it represents a single star with 
the mass of the sink particle. We determine the 
luminosity of sink part i cles u sing the calculations of 
IWuchterl fc Tscharnuterl (120031). s pecifically their Table 
3. IWuchterl fc Tscharnuterl (|2003l ) model the very early 
stages of stellar evolution. In order to calculate the lu- 
minosity of their stars, they include accretion and con- 
traction luminosity as well as stellar evolution models 
to calculate the luminosity from deuterium burning. To 
calculate the accretion luminosity, they compute the ac- 
cretion rate from fluid equations. This differs from other 
pre- main seque nce models which either ignore the mass 
accretion rate (jD'Antona fc Mazzitellil 11994 ISiess et all 
120001) or fix the mass accr et ion rate in their mod- 
els (IZinnecker fc Yorkejl2007fh IWuchterl fc Tscharnuterl 
( 2003) define a mass of the optically thick region in their 
models, M T . We assume that our sink particle mass is 
equivalent to their M T . We assume that the value of 
M, ., defined as the mass o f the optically thick region 
by IWuchterl fc Tscharnuterl (|2003f) . is equivalent to our 
sink particle mass. Therefore, given an M T and mass 
accretion rate, we can use their models to determine the 
time-averaged luminosity due to accretion and deuterium 
burning. 

In this work, we assume that objects with M < 
O.OIA/q do not have a high enough luminosity to af- 
fect the temperature of the surrounding gas; therefore, 
for these objects we set L = 0. For sink masses of 
0.01 < M < 0.04M Q , we assume that the luminosity 
is only a function of the sink mass and is independent of 
the mass accretion rate, M. We assume the luminosity 
follows the form 

Iog(L/i ) = 3.5 log(M/M ) + 5. (2) 
This relation can be deriv ed from Figure 4 of 
IWuchterl fc Tscharnuterl (|2003f) . For objects with M > 
0.04M Q , the luminosity is a function of M and M. We 
determine M at a given time by calculating the total 
mass accreted by the source over the previous 5,000 
years. 

For objects with masses greater than 2M , we need to 
calculate th e luminosity even though it is not given in 
the work of IWuchterl fc Tscharnuterl (|2003l ). We calcu- 
late the luminosity by assuming that the star is contract- 
ing on a Kelvin-Helmholtz timescale. The luminosity of 
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TABLE 2 
Luminosity Calculation 
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Derived 




L', M', M' 
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^acc 
^acc 
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= V — V 
= (M')°'^' S ' 
= GM'M'/L' m 


At 


M, M 


^M.S. 
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Lace 

L 


= M i:l 

— (R'acc - r ms )/*kh 
= i?4 cc - dR/dt x At 

= GMM/Rajx 

= Lyl.Q. + Lace 



a young star depends on the combination of luminosity 
from nuclear burning, contraction, and accretion. Mod- 
els of pre-main seque nce evolution of interm ediate- and 
high-mass stars (e.g., iPalla fc Stahler|[l993D show that, 
without accretion, young massive stars evolve at constant 
luminosity before they reach the main sequence. There- 
fore we assume that the luminosity (ignoring accretion) 
of a massive star is equal to its main sequence luminosity. 
We define the total luminosity of a star with M > 2Mq 
as a combination of its main sequence luminosity and its 
accretion luminosity, L — Lm.s. + L acc 

To calculate the main sequence lu minosity we use the 
relatio n L/L Q ~ (M/M ) 3 - 7 from [Zinnecker & Yorkd 
(|200l . We summarize our calculation in Table [H In 
order to calculate the accretion luminosity we use the 
properties of the star from the previous timestep. We 
represent these values as primed symbols, i.e., L 1 , L' acc , 



L' M s , M', M', R' acc , and R' M s , total luminosity, accre- 
tion luminosity, main sequence luminosity, mass, mass 
accretion rate, accretion radius (the actual radius of the 
object), and main sequence radius (the radius of the star 
assuming it is not accreting and has reached the main se- 
quence), respectively, all in solar units. At the timestep 
when M becomes greater than 2M Q , the known values 
are L', M' , and M'. The values that we derive are 
L acc: l 'm.s.> Kco and r 'm.s.- The main sequence lu- 
minosity is L' M s /Lq = (M'/Mq) 3 - 7 and the accretion 
luminosity is L' acc = L' — L' M s . The main sequence ra- 
dius c an be calculated using R [ f g /R P) — (M'/Mq) - 5 
(from [Zinnecker fc~Y orkc 2007). The accretion radius is 

R acc = GM'M' /L acc . These values are used to calculate 
the contraction timescale later. 

We now have all of the parameters needed to cal- 
culate L at the current timestep (un-primed symbols 
represent values at the current timestep). Therefore, 
L = Lm.s. + L acc . Lm.s. can be calculated from the 
mass. However, L acc is more difficult to calculate. We 
attempt to include the contraction luminosity through 
the contraction timescale when we calculate the accre- 
tion luminosity; for a known main sequence and accre- 
tion luminosity, we assume that a star contracts on the 
Kelvin-Helmholtz timescale, 



3M' 2 G 



(3) 



Then we can calculate the rate of change from the 
previous accretion radius to the main sequence radius, 
dR/dt = (R' acc — R'ms )/*KH- Now that we have the 



rate of change of the radius, we can calculate the new 
accretion radius, i? aC c = -R^cc — dR/dt x At, where At 
is the change in time between the previous and current 
timestep in the simulation. Then we can derive the new 
accretion luminosity, L acc = GM M / R acc . Combining 
the new accretion luminosity and new main sequence lu- 
minosity gives us the new total luminosity. As we dis- 
cuss later in §4.41 we found that this approach leads to 
a smooth transition from M < 2Mq to M > 2M , with 
no discontinuities in the luminosity. 

2.4. Dust and Gas Temperature 

The temperature of each SPH particle is initially set 
to 5K. When the first sink particle forms and has a 
non-zero luminosity, then we begin calculating the dust 
temperature. In order to calculate the dust tempera- 
tu re in our sim u lation we use the calculation described 
in lUrban et alj (|2009[ ). This method uses a spherical 
continuum radiative transfer code to calculate the dust 
temperature around individual sources with a given lumi- 
nosity and surrounding, envelope density profile. Using 
the temperature profiles around individual stars to cal- 
culate the flux at every point in the simulation, the dust 
temperature can be calculated. The gas temperature is 
then set equal to the dust temperature. 

Some differences between the method described in 
lUrban et alJ ()2009f ) and the method we use in this pa- 
per are discussed. Instead of a polynomial interpolation 
method, we use weighted-linear interpolation to decrease 
the computational time spent on t his step. We conduc ted 
tests similar to those performed in lUrban et al.l (|2009l ) to 
test our new interpolation method. We do not find signif- 
icant differences between the two interpolation methods. 

We assume that the outer radius of our dust-gas en- 
velope is O.lpc for all of the sink density profiles. This 
is because the density profile derived using the method 
described in §2.21 alw ays derived a outer radius less than 
0.1 pc. We found in lUrban et"aT1 (|2009t ) that the outer 
radius does not have a significant effect on the dust tem- 
perature. Therefore, we are not concerned that this as- 
sumption will strongly affect our dust temperature cal- 
culation. 

We have also needed to extrapolate in a few cases when 
the parameters of our sink particles wer e outsi de of the 
parameter space studied in lUrban et al.1 (|2009h . In §4.21 
and ij4.4[ we discuss the sink properties that determine 
the dust temperature - the envelope's density profile and 
the sink particle's luminosity. Based on the values of 
luminosity (i.e. L < 10 6 Lq), it is unlikely that extrap- 
olation occurred due to luminosities that were outside 
the parameter space. However, this is not the case for 
the density profile. For a > 2.0 and log n > 6, ex- 
trapolation was more frequent. But this extrapolation 
only occurred in the range of a between 2 and 2.5 and 
for valu es of log n n betwe en 6 and 8. Figures 12 and 
13 from lUrban et all (|2009l ) show that for this range of a 
and log n the behavior of K and /3 (the parameters that 
determine the dust temperature) is somewhat regular for 
low luminosities. For higher luminosities (L > 10 4 Lq), 
the behavior is not as regular. This may introduce some 
uncertainty in our temperatures around high luminosity 
sources. 

In some cases we were unable to derive a density profile 
because the sink particle was surrounded by too few gas 
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particles. When this occurred, we assumed a = and did 
not calculate a density profile for that sink. Therefore, 
at that timestep, the sinks with no density profiles were 
not used to calculate the dust temperature in the simu- 
lation. In some cases, the value of a was found to be less 
than 0.5. When this occurred, we simply assumed that 
the value of a was 0.5 and used that value to calculate 
the value of n . This assumption must be made and is 
due to our general assumption of spherically distributed 
material around sink particles. This approximation does 
not occur very often in our simulation as can be seen 
later in Figures [5] and [6l 

3. THE SIMULATIONS 

We perform four different simulations. They incorpo- 
rate all of the same physical processes (i.e., hydrodynam- 
ics and gravity) and are identical in all aspects excepting 
the following two. (1) Two of the simulations assume 
an isothermal equation of state at a fixed temperature of 
5K throughout the entire calculation (T gas = 5K). The 
other two simulations include the effect of dust heating 
discussed in §2.41 (T gas = Id us t)- For each method of 
calculating the temperature, we simulate a small and a 
large region by changing the level of particle splitting. (2) 
Two simulations allow one generation of particle splitting 
(N gen = 1) and two simulations allow two generations of 
particle splitting (N gcn = 2). The size and mass of these 
simulations depend on the number of generations of par- 
ticles splitting and the values are presented in Table [TJ 
(The values for a simulation with no particle splitting, 
Ngcn = 0, are shown only for reference.) A higher level 
of particle splitting allows larger and more massive re- 
gions to be simulated. 

Our SPH algorithm simulates the evolution of a region 
deep within a molecular cloud as a cubic volume with 
periodic boundary conditions, containing initially 64 3 , 
or 262, 144 particles. The volume of the box is L^ ox ; see 
Tableland §23] for details on ibox- Particle splitting, 
sink formation, and accretion onto sinks are included. 
For details of the SPH algorithm, we refer the reader to 
MES06 and §3 

3.1. Initial Conditions and Simulation Parameters 

Our method for generating initial conditions is de- 
scribed in detail in MES06. We start with a uniform den- 
sity distribution, with no overall density gradient. Onto 
this we superpose a small density perturbation which 
is described by a Gaussian random field with a density 
power spectrum P(k) oc fc~ 2 , where k is the wavenum- 
ber. In terms of the simulation setup, these initial con- 
ditions are achieved by arranging the SPH particles on a 
64 x 64 x 64 cubic grid, and slightly displacing them to 
reproduce the desired power spectrum. 

To generate initial conditions, we first need to fix the 
mass resolution of the algor ithm. In their SPH simula- 
tions, |BatiI&^onne3 (|2005[ ) use a barotropic equation of 
state, which is isothermal at densities p < 10 _13 gcm~ 3 
(or n « 2.6 x 10 10 cm -3 ), and adiabatic at densities 
p > 10~ 13 gcm -3 . In this case, there is a minimum Jeans 
mass (Mj) m ; n = 0.0011 M©, corresponding to the den- 
sity p = 10~ 13 gcm~ 3 . Since we focus on the formation 
of a star cluster and not the details of individual star 
formation, it is unnecessary for us to resolve such high 
densities and low masses. Instead, we set the mass res- 



olution limit of our algorithm at M = 0.008 M© w 10% 
of the minimum mass for hydrogen-burning. This is not 
much of a limitation, since the effect of dust becomes im- 
portant only at much larger masses/higher luminosities. 
A drawback of this limited resolution is that we will not 
be able to assess the effect of heating by massive stars on 
the subsequent formation of objects with masses below 
0.008 M©. However, this lower resolution allows us to 
simulate a large region which contains enough mass to 
potentially form several high- mass stars. 

We consider a cloud with initial density p = 4.75 x 
10 _20 gcm -3 , or n = 1.22 x 10 4 cm~ 3 assuming p = 2.33 
(p = /inmn), and temperature T = 5K (see Table \5§ . 
Our choice of an initial tempe rature of 5K, as o pposed to 
10K used in other simulations (jBate et alj2003t MES06), 
was mo tivated by th e low temperatures in the calcula- 
tions of Urban et alj (|2009fl as well as the discussion in 
ILarsonl ( 2005) and recent observ ations from lEvans et all 
(1200 Hi and ICrapsi et ail (|2007j) . The gas will remain 
isothermal at 5K, until it is heated by dust. 

Our algorithm will turn dense gas clumps of mass M = 
0.008 M© into sinks. To justify this, these objects must 
have a mass equal to the Jeans mass. For a gas with 
polytropic constant 7 = 5/3, the Jeans mass is given by 

(|Tohlind [T982K Since the gas is isothermal, we can set 
T = 5K and Mj = 0.008 M© in equation @, and solve 
for the density. We get p c — 2.822 x 10~ 16 g cm~ 3 , or 
n c = 7.252 x 10 7 cm~ 3 . This is the threshold density 
at which sinks will be created (see Table [3]) . This will 
happen after the gas contracts from the initial density, p, 
by a factor of p c jp — 5, 942. This factor is smaller than 
the value of 40,000 used in MES06, but still provides a 
wide dynamical range in density. 

It takes a minimum number of particles to properly re- 
solve a Jeans mass. The precise value depends on the par- 
ticular implementation of SPH. In MES06, we assumed a 
typical value of 100 particles. In this paper, we also split 
particles when the Jeans mass drops below 100 particles, 
but we require that the clumps that turn into sinks con- 
tain 200 particles. The reason is that clumps with 100 
particles often fail the criteria for sink creation because 
their rotation or internal motions make them unbound. 
Therefore, the mass of an SPH particle is 

mp „ = lS = 4xl0 - 5 „ s . (5) 



200 



Our simulations start with 64 3 particles. Therefore, the 
total mass of the system is M tot = 64 3 m par t = 10.49 M©, 
and the box size is ibox = (-Mtot/p) = 0.246 pc. By 
setting p = p and T = 5K in equation (01, we get an 
initial Jeans mass Mj : [ nit — 0.617 M©. Hence the system 
starts up with Nj = 17 Jeans masses, compared to 500 
in MES06. However, these numbers assume no particle 
splitting. If we allow particles to split, each splitting gen- 
eration will increase the effective number of particles by 
a factor of 8. The densities at which particle splitting 
occurs are listed in Table [31 Since the mass of the final 
generation of particles is fixed by equation |(5J), each split- 
ting generation will increase M to t and Nj by a factor of 
8 and Lbox by a factor of 2. With a fixed density resolu- 
tion and a fixed minimum mass resolution, we effectively 
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TABLE 3 

Density Parameters 



Density g cm s cm 

Initial average density (pi) 4.75 X 10" 20 1.22 X 10 4 

Sink creation density 2.82 X 10~ 16 7.25 X 10 7 

Splitting Density for N gon = 1 371 X p.; 

1st and 2nd Splitting Density for N ge n = 2 5.8 X pi 371 X pi 



increase the size and mass of our simulated region with 
each new level of particle splitting. To get systems with 
reasonable sizes and masses, we allow for 1 and 2 gener- 
ations of particle splitting, N gon = 1 and 2. 

As mentioned before, our density resolution is not 
as high as that of MES06. If we compare to Case 4 
listed in Table 4 of MES06, we find that they start 
with a higher initial density and model a smaller region 
(Lbox = 0.38pc) with a larger mass (M to tai = 320M©) 
than our N gen = 1 model. However, our N gen = 2 model 
is larger in volume and more massive, but still starts with 
a smaller initial density. 

3.2. Sink Particles 

As mentioned in the previous section, sink particles are 
formed from 200 gas particles which exceed the threshold 
density for sink creation (p c = 2.822 x 10~ 16 g cm -3 ). 
The second criterion for sink creation is that the gas 
particles must also be Jeans unstable, meaning that they 
are collapsing and have formed a gravitationally bound 
system at the time that the sink is created. The exact 
prescription for creating a sink particle is described be- 
low. 

We define an accretion radius r acc , such that when a 
particle reaches the threshold density p c , a sphere of ra- 
dius r acc , centered on the new location of the sink parti- 
cle, will contain 200 gravitationally bound particles. In 
the isothermal simulations, this corresponds to a Jeans 
mass. These 200 particles are then removed and replaced 
by a sink particle with the same total mass and center-of- 
mass position and velocity as the 200 particles. In order 
to determine the value of the accretion radius, we have 
run test cases in the isothermal simulation in which we 
vary the value of the accretion radius and allow a few 
sink particles to form. We choose the value of the ac- 
cretion radius that ensures that the number of particles 
that are used to create a sink particle is approximately 
200. For our simulations the accretion radius is set at 
~ 150AU. 

For the runs that include dust heating, the tempera- 
ture will vary throughout the simulation so there will be 
no fixed value of the Jeans mass. In this case, we still 
use the accretion radius from the isothermal simulation. 
The first criterion of needing ~ 200 particles will still be 
met if we use the isothermal accretion radius. However, 
the second criterion requiring that the particles must be 
Jeans unstable before a sink forms will tend to delay sink 
formation compared to the isothermal simulation. This is 
because the Jeans mass in a simulation with dust heating 
will typically be higher due to the increased temperature. 
Hence, it will typically take more than a collection of 200 
Jeans-unstable gas particles to have a total mass equal 
to the Jeans mass; when a particle reaches the thresh- 
old density, sink formation will be delayed until enough 



particles (or mass) are within the accretion radius to in- 
crease the enclosed mass to the local Jeans mass defined 
by the local temperature. 

Once sink particles are formed, they have the ability to 
grow by accreting gas particles. Whenever a gas particle 
enters the accretion radius of a sink, that particle will 
be accreted, provided that it is gravitationally bound to 
the sink. Typically the sink particles will have time to 
accrete mass before their luminosity is large enough to 
heat the surrounding dust and gas. 

3.3. Timescales 

Since our initial density is p = 4.75 x 10 _20 gcm~ 3 , or 
n = 1.22 x 10 4 cm~ 3 , the initial free-fall time, defined as 

i ff = y/3n/32Gp, (6) 

is t s = 9.64 x 10 12 s = 3.06 x 10 5 yr. This is true for all 
simulations because the initial density is the same for all 
cases. 

Since we do not include ionizing radiation in our sim- 
ulations, they are no longer realistic when very massive 
stars form. Therefore, we stop our simulations when the 
most massive sink particle in each reaches a mass of 20.8 
Mo, which is approximately the mass of an 09.5 star. 
iKetol (|2003f ) shows that for earlier spectral types, the 
HII region surrounding a star (or group of stars) could 
be greater than the size of our sink particle radius or 
would be expanding. It is no longer realistic to ignore 
the effects of ionization. Therefore we stop our simu- 
lation when we have reached this limit. This occurs at 
2.5tff for the N gon = 2 simulation with dust heating. We 
stop our N gcn = 2 isothermal simulation at the same time 
so that we can compare the two N gcn = 2 simulations. 
At t = 4.5tff, in our N gcn = 1 simulation, no sinks have 
yet reached 20.8 M . However, over 80% of the mass is 
in sinks. Therefore, we stop the simulation at this point. 

3.4. Gas Temperature 

Our version of SPH does not include the standard en- 
ergy equation, with p dV work, viscous heating, and ra- 
diative cooling. Instead, the temperature of each particle 
is set at each timestep, depending on the local conditions. 
In the isothermal runs, the temperature is kept fixed, 
while in the runs with dust heating, the gas temperature 
is set to the local temperature of the dust (see £|2.4p . 
We run two isothermal simulations, and two simulations 
with dust heating, both with one and two generations 
of particle splitting, for a total of four simulations. In 
Table |4l we list the four different simulations, and give 
the final number of sinks and the final mass in sinks. 

4. RESULTS 

Figures [1] and [2] show how the mass in the simulations 
is distributed between the gas and the sink particles as 
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TABLE 4 
Summary of Simulations 



Agcn Tgas Final Time Final # of Sinks Max. Sink Mass Total mass in Sinks SFR ff 
(iff) (Mq) % 



1 5K 

2 5K 

2 Td ust 



4.5 
4.5 
2.5 
2.5 



518 
20 
3429 
74 



0.46 
12.4 
0.50 
20.8 



98% 
87% 
60% 
50% 



0.22 
0.19 
0.24 
0.20 



1 1 1 1 1 1 1 1 1 1 1 . 




- T x= 5K \ 




Cbs Fraction \ 




Sink Fraction 








J 






■ 

o.b 

m 0.6 

f 

0.4 
0.3 



■ 


- VT^ "\ 




O.B 






n 0.6 


Gas Fraction 




r 


Sink Fraction 




O.'l 


N - ta 




0.2 













a ■ — ' — ' — ■ — ' — ' — ' — ' — ' — ^ — I — 1 — ■ — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — ^ a 

1 £ tlm. (l„) 3 4 



Fig. 1. — Mass fraction evolution for simulations with N gcn = 1. 
Panels show fraction of mass in gas and sinks as a function of 
free-fall time. Top panel shows results for the isothermal simula- 
tion. Bottom panel shows results for simulation with dust heating. 
Thick line shows the mass fraction in sinks. Thin lines show the 
mass fraction in gas as a stacked histogram for different genera- 
tions of particles. Lower thin line represents particles which have 
not undergone particle-splitting. Higher thin line represents frac- 
tion of particles which have split once. The two gas lines are barely 
indistinguishable for the simulation with dust heating. Dashed line 
shows the evolution of the number of sinks. Vertical line at 2.4tg 
is shown for reference. 

a function of time. These figures show that as the mass 
in gas particles decreases, the mass in sink particles in- 
creases, as expected. It is also interesting to note that 
the transition from a gas-dominated to sink-dominated 
simulation (point at which the sink and gas lines cross in 
Figs. Q]and[2]) occurs at different times for the different 
simulations . Since dust heating increases the average 
temperature of the simulation box (see ij4.3l below) and 
thus inhibits the formation of sink particles, the tran- 
sition from gas-dominated to sink-dominated occurs at 
a later time in simulations with dust heating, both for 
Ngcn = 1 an d N gen = 2. Indeed, for the case N gen = 2 
with dust heating, the mass in sinks has not yet reached 
50% by the end of the simulation, at t = 2.5iff. 

The increase in temperature that inhibits the forma- 
tion of sink particles also affects particle splitting. In 
the simulations with dust heating, the higher tempera- 
ture leads to a higher local Jeans mass. Therefore, the 
condition which triggers particle splitting, i.e. violating 
the Jeans criterion, is rarely met at the highest level of 
particle splitting for the simulations with dust heating, 
unlike the case for the isothermal simulations. This can 
be seen in the lines representing "Gas Fraction" in Fig- 
ures [T] and [2l as well as the lack of blue and green dots 
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Fig. 2. — Mass fraction evolution for simulations with -/V gcn = 2. 
See Figure[T]for details of plot. The line representing particles that 
have split twice is barely indistinguishable from the line represent- 
ing particles that have split once in the bottom figure. 

in the right-hand panels of Figures [3] and [H respectively. 

Another interesting feature of Figures [1] and [2] is the 
large difference in the number of sink particles formed in 
the simulations. This can also be seen in Tables [5] and [51 
The isothermal simulations form more than an order of 
magnitude more sink particles than the equivalent simu- 
lations with dust heating. However, as seen in Table |4j 
the amount of material in sinks is comparable for simu- 
lations with similar sizes, i.e., same values of N gen . We 
discuss these two features next. 

The difference in number of sinks formed is strongly 
affected by dust heating. As the sink particles heat the 
simulation volume in the simulation with dust heating, 
the growth of structure and sink particle formation is 
inhibited. The heated gas prevents the fragmentation 
which occurs unhindered in the isothermal simulation. 
Since the percentage of the total mass in sinks is com- 
parable for simulations with the same value of N gen and 
there are fewer sinks in the simulations with dust heat- 
ing, the sinks in the simulations with dust heating are on 
average much more massive than the sinks in the isother- 
mal simulations. Gas that is prevented from forming new 
sinks instead accretes onto existing ones. This will be 
discussed in §4.41 below. 

The percentage of material in sink particles at the end 
of the simulations, or the star formation efficiency, is 
affected slightly by the energetics algorithm used, i.e. 
isothermal or with dust heating. As seen in Table |4] 
in the N gen = 2 calculation, 60% of the mass is in sinks 
after 2.5 free-fall times for the isothermal calculation ver- 
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TABLE 5 
Simulation Summary: N k , 




TABLE 6 
Simulation Summary: N e , 



time 


T gas = 5K 








(%) 


Max. Mass (M Q ) 


# of Sinks 


Max. Mass (M ) 


# of Sinks 


1.6 


0.04 


14 


0.066 


3 


1.8 


0.11 


899 


2.11 


28 


2.0 


0.20 


2452 


7.99 


53 


2.2 


0.34 


3033 


12.40 


70 


2.4 


0.46 


3371 


17.99 


71 


2.5 


0.50 


3429 


20.8 


74 



sus 50% for the simulation with dust heating. Another 
measure of the star formation effici ency (or the speed of 
star fo rmation) has been defined bv lKrumholz fc McKed 
(2005). SFRff is the fractio n of material converte d 
into stars per fr e e-fall time (|Krumholz fc Tanl I2007D . 
IKrumholz fc~Tanl (pOOl argue that SFRff ranges over 
0.013 - 0.03 for the Milky Way Galaxy based on the inter- 
pretation of vari ous observations, whi c h is in agreement 
with the work of IKrumholz fc McKeel (|2005f h Studies of 
nearby clouds forming low-mas s stars found an average 
value of 0.04 (jEvans et al.ll2009fh 

We give the values of SFRff for our simulations in Ta- 
ble 2] All of our results are abo ut an order of magnitude 
higher than the prediction of IKrumholz fc TanT (120071) 
and a factor of 5 higher than the results of I Evans et al.l 
( 2009). The values of SFRff for the simulations with dust 
heating are slightly lower than the values for the isother- 
mal simulations. We believe the higher values of SFRff 
in our simulations probably arise b ecause we are not in- 
cludin g turbulence. In the theory of lKrumholz fc McKeel 
(2005), the main regulating agent of star formation is 
turbulence. Since we do not include this physical effect, 
which would slow star formation in our simulation, it is 
not surprising that we find higher values of SFRff . 

4.1. Sink Particle and Gas Mass Distribution 

Figures [3] and Q] show the distribution of sink and gas 
particles in our simulation box. They are shown at t = 
2.4 tg because at this time a substantial fraction of the 
total number of final sink particles have formed and there 
is also still a significant amount of gas remaining, as seen 
in Figures [T] and [51 Therefore it is an appropriate time 
to compare how different equations of state affect sink 
particle formation. 

The most noticeable difference between the isothermal 
simulation and the simulation with dust heating is the 
sparseness of sink particles in the simulation with dust 
heating. Less particle splitting is also occurring. When 



the first sink particles begin to heat the environment, the 
Jeans mass increases and fragmentation is halted. An- 
other consequence of dust heating is the lack of definition 
in the filaments for the simulations with dust heating 
compared to the isothermal simulation (seen clearly in 
Fig. [31) ■ The hotter temperatures prevent the filaments 
from collapsing toward the high density central region. 
In Figure [4[ there are more filaments in the isothermal 
simulation (left-hand panel) because fragmentation pro- 
ceeds unhindered by any increase in temperature, which 
occurs in the simulation with dust heating (right-hand 
panel). 

4.2. Density Profile Evolution 

At each time step in the simulations with dust heating, 
we calculate the density profile around each sink particle. 
As the sink particles move around in the simulation box 
and material accretes, their surrounding density fields 
change. The two parameters that define the density pro- 
file, n Q and a are discussed in 32.21 Figures [5] and [6] show 
the evolution of the density profile for the sink particles 
formed in the simulations. Figures [7] and [H] show the aver- 
age values and dispersions of a and n for the individual 
sink particles. 

We calculate the average values of a and n a using vari- 
ous methods. The data are sampled from the simulation 
every 0.05£ff. Of the points plotted in Figures [5] and [6l 
we calculate the average value of (1) all of the points, 
(2) only points with a ^ 0, and (3) only points with 
M > lO~ 8 M0yr~ 1 . We summarize all of the results in 
Table [7] The spread is very large when we consider all 
points. Ignoring sinks that are not likely to be accret- 
ing, i.e., those with a = and with low accretion rates, 
we find smaller dispersion in the results. For case (3), 
we find (a) = 1.7 ± 0.3, (log(n /cm- 3 )) = 6.2 ± 0.3 
for simulations with N gon = 1 and (a) = 1.7 ± 0.4, 
(log(n /cm~ 3 )} = 6.5±0.3 for simulations with N gcn = 2. 

We can compare our average density profile values 
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Fig. 3. — XY position plot of sink and gas particles for simulations with N gGn = 1 at 2.4 tg. Left plot shows result from isothermal 
simulation. Right plot includes dust heating. Black and blue dots indicate gas particles. Blue dots are gas particles which have undergone 
one particle splitting. Red dots are sinks. Scale is 0.492 pc. X 0.492 pc. 




X X 

Fig. 4. — XY position plot of sink and gas particles for for simulations with N gGn = 2 at 2.4 iff. Left plot shows result from isothermal 
simulation. Right plot includes dust heating. Black, blue, and green dots indicate gas particles. Blue dots are gas particles which have 
undergone one particle splitting. Green dots are particles which have split twice. Red dots are sinks. Scale is 0.984 pc X 0.984 pc. 



to those derived observationally. Class and Class 
I cores, which are representative of the earliest stages 
of isolated low-mass s tar fo rmation, have been studied 
bv lShirlev etHl rt20f)2h an d lYoung et all (I2003T ). respec- 
tivelv. IShirlev et all (|200l find (a) = 1.63 ± 0.33 and a 
typical value ofa=1.8±0.1if the y ignore two s ources 
with aspherical emission contours. lYoung et all (I2003T ) 
find (a) = 1.6 ± 0.4. There is excellent agreement be- 
tween the average values of a derived for the density 
profiles around sinks in our simulation and the observed 
values of a in isolated low-mass star-forming cores. The 
values of n a derived from the tables in the papers give 
(log(n /cm- 3 )} = 6.1 ± 0. 2 (IShirlev et all l2002f) and 
_(log(n /cm- 3 )} = 5.4 ± 0.5 (jYoung et al.ll2003l ). There 
is some agreement with our average values; however this 
comparison may be affected more strongly by other fac- 



tors such as the masses of the individual cores and the 
observed intensity to density conversion. 

Another interesting feature of the density profile pa- 
rameters is the relationship between the dispersion of a 
and n a and the final mass (shown in Figure [5]) . We find 
that objects with high dispersions tend to be the lower 
mass objects (M < 5M Q ) in our system. Conversely, 
objects with the highest mass are more likely to have 
low dispersions. (This trend does not appear to apply 
for sinks formed after 2tg; however, these sinks may not 
have had enough time to accrete 5M©.) 

4.3. Temperature and Density Evolution 

Several groups that model clustered star formation 
assume a s implifi ed equation of state. For example, 
iBate et all (|2003l ) assume that the gas is isothermal up 
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TABLE 7 
Density Profile 



Case 


AT gcn = 1 




Ar gc „ = 2 






(a) 


(log(ra /cm 3 )} 


(a) 


{log(n /cm 3 )} 


All points 


0.8±0.9 


2.9 ±3.1 


1.6±0.6 


5.8±2.0 


Points with a > 


1.7±0.4 


6.1 ±0.3 


1.8±0.4 


6.5±0.3 


Points with M > 10~ 8 M Q /year 


1.7±0.3 


6.2 ±0.3 


1.7±0.4 


6.5±0.3 
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Fig. 5. — Evolution of the density profile parameters for simula- 
tion with N gcn = 1, including dust heating. Each color/line type 
represents the evolution of a different sink. Data are sampled every 
0.05t ff . 
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Fig. 6. — Evolution of the density profile parameters for simula- 
tion with N gcn = 2, including dust heating. Each color/line type 
represents the evolution of a different sink. Data are sampled every 
0.05%. 

to a density of 10 _13 g/cm 3 (or n = 2.6 x 10 10 cm~ 3 ). 
However, this work ignores t he eff ect of stellar heating 
of the dust and gas. iLarsonl (|2005f ) gives an equation of 
state which represents the state of the gas before stars 
are born and have significant luminosity. We show the 
form of this equation as a solid line in Figures [9] and 



Fig. 7. — Density profile parameters for simulation with N geI1 = 
1, including dust heating. Values of density profile parameters, a 
and n , are shown for the sinks. The time-averaged values of a and 
n over all sinks are given in the figure. Error bars shown with a 
solid line indicate the standard deviation for each individual sink. 
Error bars shown with a dotted line indicate the minimum and 
maximum value of a and n a . 
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Fig. 8. — Density profile parameters for simulation with N sen = 
2, including dust heating. Similar to Figure[7] The points that are 
red are sink particles with final masses greater than 5Mq . 

[TU1 These figures also show the temperature and density 
of the gas particles at various times in our simulations. 
The fingers that are seen extending to the right in these 
plots correspond to the gas particles that are close to a 
sink particle and therefore have their dust temperature 
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TABLE 8 

Temperature (K) Statistics for N scn = 1 

time (iff) mean sigma median lower upper 

quartile quartile 



2.0 


5.6 


2.1 


5.0 


5.0 


5.0 


2.5 


24 


3.1 


22 


20 


26 


3.0 


33 


13 


29 


26 


34 


3.5 


36 


11 


32 


30 


38 


4.0 


37 


15 


32 


30 


38 


4.5 


39 


20 


31 


28 


10 



TABLE 9 

Temperature (K) Statistics for N gan = 2 

time (iff) mean sigma median lower upper 

quartile quartile 



1.6 


2.4 


1.9 


5.0 


5.0 


5.0 


1.8 


15 


6.1 


12 


11 


16 


2.0 


24 


9.4 


21 


19 


25 


2.2 


34 


15 


29 


26 


34 


2.4 


47 


29 


38 


33 


19 
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Fig. 10. — Temperature versus density relation of SPH particles 
as a function of time for simulation with iV gcn = 2. The statistics of 
the temperature are listed in Ta ble[9] Solid line shows the equation 
of state given by ILarso 3 ll2005h . 
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Fig. 9. — Temperature versus density relation of SPH particles as 
a function of time for simulation with N gcn = 1. The statistics of 
the temperature are listed in Ta ble[8] Solid line shows the equation 
of state given by Larson (2005). 

determined by the luminosity of a single sink particle. It 
is clear from Figures [5] and [TU] that a simple equation of 
state that describes the behavior of the temperature as 
a function only of density is insufficient once stars begin 
to form. When this happens, heating depends not only 
on the local density, but also on non-local effects, such as 
the distance to a luminous sink and the star formation 
history. 

We give temperature statistics for the two simulations 
in Tables [8] and [9] For both simulations with dust heat- 
ing, we find that there is a general trend for denser gas 
to be hotter. Another interesting point is the behavior 
of the lower quartile which gives the median value of the 
lower 25% of the distribution. This value increases as 
a function of time indicating that overall the entire re- 



gion is getting hotter. As the sources turn on they slowly 
heat all of the material in the entire volume. This "global 
warming" is what changes the time evolution and the re- 
sulting mass function compared to the isothermal case 
(see SH) . 

4.4. Luminosity and Mass Evolution 

Figures [11] and show the luminosity and mass ac- 
cretion rate as a function of mass for the two simulations 
which include dust heating. (Recall that the mass accre- 
tion rate is averaged over a time interval of 5000 years, see 
§2.30 The luminosity as a function of mass follows the 
same trend in both simulations. They both also follow 
the main sequence mass-luminosity relation, L oc M , 
at high masses with additional luminosity due to accre- 
tion. Sink particles that suddenly stop accreting can be 
seen by the lines that drop sharply from the trend of the 
rest of the sink particles in the plot of mass accretion rate 
versus mass. When a sink particle has a low luminosity 
it can be seen in the mass accretion plot to be in a phase 
of low accretion. (The same holds for high-luminosity 
sink particles and high accretion rates.) Sink particles 
shown by the red and green dashed lines in Figure [TT] 
at ~ O.IMq demonstrate this effect for low- mass sink 
particles. 

These figures also sho w that our method of calculat - 
ing the luminosity using IWuchterl fc Tscharnuterl (|2003l ) 
and the method described in §2.31 are compatible with 
one another; the luminosity varies smoothly at 2Mq (the 
transition mass between the two methods). There is a 
slight increase in the slope at 2M . This is probably due 
to the assumption that the stars with M > 2Mq have at 
least a main sequence luminosity as discussed in §2.31 

In Figure I13[ we show the relation between luminosity 
and mass of the individual sinks in separate panels. This 
figure demonstrates that high-mass sink particles have a 
larger fraction of their luminosity due to their main se- 
quence luminosity as they become more and more mas- 
sive. 

There are interesting differences between the two sim- 
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Fig. 11. — Luminosity and mass accretion rate as a function 
of mass for simulation with JVg en = 1 and dust heating. Values 
of mass, mass accretion rate, and luminosity are plotted for all 
sinks created in the simulation with N gen = 1 and dust heating. 
Different colors and line types correspond to different sinks. Solid 
black line in the top panel near the top right corresponds to a slope 
of 3.7. 
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Fig. 12. — Luminosity and mass accretion rate as a function of 
mass for simulation with -/V gcn = 2 and dust heating. See Figure 
1 1 1 1 for details. 

ulations with dust heating. For the simulation with 
N gcn — 1 (Fig. Hip , the maximum luminosity, mass, and 
accretion rate are all lower than the maximum values in 
the simulation with N gcn = 2 (Fig. [T2"|) . The fact that 
the -/Vg Cn = 2 simulation has larger values of luminosity, 
mass, and accretion rate is likely to be due to the larger 
total mass inside the simulation volume (Table [TJ. The 
sink particles have more available mass to accrete, which 
leads to higher sink masses which in turn affects their 
mass accretion rate and luminosity. Another way to in- 
terpret this result is that the higher mass in the AT gcn = 2 
simulation leads to more of the rare, large fluctuations 
in the density field which allow bigger objects to form. 

Figures [T4l and [T5l show the mass, luminosity, and mass 
accretion rate as a function of time for the two simula- 



FlG. 13. — Luminosity evolution as a function of final mass 
for the dust heating simulation with N gcn = 2. Top row shows 
the luminosity evolution of sink particles with M > IOMq. Sec- 
ond row shows evolution of sink particles with masses between 1 
and 10 Mq. Bottom three rows show evolution for sink particles 
with M < IMq. Numbers within boxes indicate order of sink 
formation. Solid line in d icates the luminosity contribution from 
Wuchtcrl & Tscharnutcr (2003) luminosity calculation and accre- 
tion luminosity, i.e. luminosity from L = M" 3,7 is not included. 
Dotted line indicates the total luminosity evolution. When there 
is no dotted line, there is a negligible contribution from the main 
sequence luminosity 

tions with dust heating. If we only compare the time pe- 
riod in our iVg en = 1 simulation which corresponds to the 
entire runtime of our N gcn = 2 simulation (i.e. — 2.5tg), 
then we find similarities and differences between the two 
simulations. As Figures 1 and 2 show, at t = 2.5ts a 
comparable percentage of the mass has been converted 
into sinks for both simulations, even though the number 
of sinks formed is different, 10 sinks for N gcn — 1 and 74 
sinks for N gea — 2. Comparing Figures [T4l and [TBI during 

— 2.5ifT, the formation time between sink particles is 
greater in the N gcn = 1 simulation. This is due to the 
larger sample volume in the N gon = 2 simulation. 

Another interesting difference in the small and large 
simulations is the time of formation of the most massive 
object. Although the most massive object in the N gcn = 

1 simulation only reaches 12.4 Mq, it is one of the first 
objects to form in the simulation. This is not the case for 
the iVgon = 2 simulation. For this simulation, the most 
massive object forms sometime after the formation of the 
10th sink particle. This can be seen in Figure 1151 The 
most massive sink particle is represented by the black 
solid line which forms at time ~ 1.72tg. 

4.5. Mass Accretion Evolution and Supersonic Accretion 

Figures [TH] and [T7] show the average mass accretion 
rates for each sink particle from the simulations with 
dust heating. We find the mass accretion rat es to be 
highly variable. This was seen previously in iKlessenl 
(|2001f ). The time-averaged mass accretion rate over all 
sink particles is M = (6.10 ± 11.47) x 1O" 6 M yr" 1 
for the simulation with N gcn = 1 and M = (2.20 ± 
3.08) x 1O~ 5 M yr" 1 for the simulation with iV gon = 2. 
These mass accretion rates are generally higher than 
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Fig. 14. — Mass, luminosity, mass accretion history of sinks 



for N E 



1. Time evolution of mass, luminosity, mass accretion 



heating. Different colors and line types correspond to different 
sinks. 
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Fig. 15. — Mass, luminosit y, m ass accretion history of first 28 
sinks for N gon = 2. See Figure ITU for details. 

those determined observationally from UV, optical, and 
IR emission excesses in classical T Tauri stars, i.e., 
M < lQ- 6 M r , yr" 1 (lHartmannll2001l lBro wn & Chandler! 
119991 I White fc Hillenbrand! l200l . However, these data 
come from objects that are more evolved than the ob- 
jects formed in our simulation and have lost most of 
their initial envelope. They are also low-mass objects 
(M < IMq). Currently, there are only a few ob- 
servat ions of mass accret ion rates for young, massive 
stars. IZapata et all (|2008f) find a mass accretion rate of 
4 — 7 x 10~ 2 M Q yr _1 for W51 North from the observed 
CN line profile. The mass of the central object, which 
could be an O star or a group of B stars, is ~ 40M©. 
IZapata et afl (|2008h also list a range of other observed 
mass accretion rates and masses. The rates derived for 
gas masses with 200-300 M Q and proto-stellar masses 
with 20-40 Mq are 10~ 4 - 10~ 2 M Q yr" 1 . These are 



Fig. 16. — Mass accretion rates in simulation with N gan = 1. 
Data was sampled at increments of ~ 0.02ig . The average value of 
mass accretion rate over all sinks is M = (6.10 ± 11.47) X 10 — 6 Mq 
yr — 1 . Error bars shown with a thick line indicate the standard 
deviation for each individual sink. Error bars shown with a thin line 
indicate the minimum and maximum value of the mass accretion 
rate. 




Sink Number 



Average value of mass accretion rate over all sinks and all time 
is M = (2.20 ± 3.08) X 10 _5 M Q yr" 1 . Thick- and thin-lined er- 
ror bars represent similar values as those discussed in Figure [TrJl s 
caption. 



higher than the high mass accretion rates seen in Fig- 
ures [16] and [T71 but the most massive object that we 
form is only 20.8Mq. 

The results shown in Figures [T|5] and [T7] are plotted 
against the Sink Number which indicates the order in 
which the sink particles were formed. This shows that 
the order of formation of the sink particles does not have 
a strong effect on the accretion rate. For the simulations 
with iVg en = 1, in which there is less mass, the accretion 
rate does appear to drop for sink particles formed later. 
However, this effect is not seen in the simulation with 
Ng en — 2, where there are still sink particles with high 
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log (Final Mass (M Q )) 

Fig. 18. — Mass accretion rates in simulation with N gen = 1 as 
a function of the final sink mass. Thick- and thin-li ned error bars 
represent similar values as those discussed in Figure IT6l s caption. 
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Fig. 20. — Sound speed (c s ) versus relative speed (v) of accreted 
particles plotted for N gen = 1 simulation with dust heating. Solid 
line shows the relation for M. = 1. 




log (Final Mass (M s )) 

Fig. 19. — Mass accretion rates in simulation with N gcn = 2 
as a function of the final sink mass. Thick- and thin-lined erro r 
bars represent similar values as those discussed in the Figure IT6l s 
caption. 

accretion rates at high Sink Number. 

In Figures [15] and [TH] we find a correlation between 
average accretion rate and the final sink particle mass 
at masses above ~ 2M Q . This suggests that high-mass 
objects are built up with large accretion rates. This 
has been suggested by others as a method of overcom- 
ing the radiation pressure wh ich could halt the forma- 
tion of massive young stars (|Yorke fc Sonnhalteil 120021 
iKrumholz etal] 120091 . This correlation was either not 
seen in other similar wor ks because they did n ot form 
objects larger than 2M (jBate fc Bonnelll 120051 ) or the 
correlation was we ak because they formed only a f ew ob- 
jects above 2M Q (|Batell2009al lOffner et all 12009 ). Had 
these earlier works started with larger initial masses, they 
most likely would have seen a stronger trend once they 
began forming more massive sink particles. 

In Figures [20] and [21] we show the sound speed ver- 




log v (km/a) 



log v (km/i) 



Fig. 21.— Similar to Figure [20] but for N gcn = 2. 

sus the speed of the sink particle relative to the accreted 
particle for the simulations with dust heating at a mo- 
ment just before the particle is accreted. Particles are 
accreting supersonically on average. The average sound 
speed of the particles increases over time due to the in- 
crease in temperature as seen in Figures H] and [TU] The 
average Mach number for accreting particles, M, can be 
calculated for our simulations. For the isothermal simu- 
lations, M = 38.4±31.5 (JVgen = 1) and M = 32.1±29.6 
(N gen — 2). The values for the simulations with dust 
heating are much smaller, M. = 11.8 ± 7.3 (N gen = 1) 
and M. — 10.8±7.4 (N scn = 2). In all cases, the accretion 
onto sink particles is supersonic; however, the velocity of 
accreting particles relative to the sound speed decreases 
when dust heating is included. 

4.6. Mass Function 

Figures l2"2T[2"5l show the mass functions of the sink par- 
ticles in all of our simulations at different times. Since 
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Fig. 22. — Mass histograms for simulation with jV gen = 1 and 
isothermal equation of state. Mass fun ction shown a t different 
tim es listed in top r ight corner of boxes. [Salpcter ( 1955) (straight) 
and Chabricr (2003) (curved) analytic mass function shown as solid 
lines. Dashed line shows initial Jeans mass in the simulation. 



1000 r 




Fig. 23. — Mass histograms for simul atio n with Ng Cn = 2 and 
isothermal equation of state. See Figure l22l for details. 

our simulations with dust heating are isothermal until 
the sink particles begin heating the surrounding environ- 
ment, the first sink particles will form at the same time 
for both the isothermal simulation and simulation with 
dust heating if the simulations have the same number of 
particle splittings (i.e., same values of iVg en ). The first 
sink particle for the simulations with iVg Cn = 1 forms 
at t ~ 2tg, but for the simulations with iVg en = 2, the 
first sink particle forms at t ~ 1.5fg . Since there is more 
mass in the simulations with N gcn = 2, there is a higher 
probability of forming a sink particle earlier. 

The isothermal simulations of MES06 all have the same 
initial mass but have different values of N gcn . Hence, an 
increase in N gen resulted in a higher resolution since the 
mass per particle was smaller. These simulations pro- 
duced log-normal distributions, with an average value 
which depended on resolution. As the resolution in- 
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FlG. 24. — Mass Histo gra ms for simulation with iV gC n = 1 and 
dust heating. See Figure [22] for details. 
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FlG. 25. — Mass histo gram s for simulation with N gcn = 2 and 
dust heating. See Figure [22] for details. 

creased, the mean shifted to lower values. In our simula- 
tions, discussed in this paper, we use a different approach 
by keeping the resolution fixed. Hence, the total mass of 
the system increases with the value of AT gen . Based on the 
results of MES06, we expected the mean value of mass in 
our simulations to be independent of N gcn . As Figures 
l22l and [23] show, this is indeed the case. We find a log- 
normal distribution with an average value of ~ 0.1 M© 
for both of our isothermal simulations. In the isothermal 
simulations, we also find that we are unable to create 
any objects with masses greater than ~ IM@. This is 
the case for many other isothermal simulations and was 
our motivation for using a dust heating algorithm. 

The mass functions for simulations which include dust 
heating (Figs. I2"4l and |2"5|) show that we are able to form 
massive stars (M > IOMq), unlike the isothermal simu- 
lations. However, for the simulation with N gcn — 1 (Fig. 
I24p , the distribution is very sparsely sampled. Adding 
more mass to the simulation leads to a more well-sampled 
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Fig. 26. — Initial versus final mass for isothermal simulation 
with iVgcn = 2. The smalles t po ssible initial sink particle mass 
is 0.008-Mq, as discussed in i|3.1l This can only increase by the 
addition of individual gas particles of mass 4 X 10~ Mq. This ex- 
plains why sink particles will only form at specific intervals when 
they form with masses close to the resolution limit of the simula- 
tion. 

distribution, even after only 2.5tg as seen in Figure [55] 
where N gcn = 2. 

In Figures [26l[27l we compare the initial and final 
masses of the sink particles. The initial Jeans mass of 
our simulations is 0.617M Q . If objects were collapsing to 
the scale of the Jeans mass without fragmentation, then 
we would expect the average mass of the simulation to be 
the initial Jeans mass. For the isothermal simulations, 
fragmentation is a process which is only halted by the 
resolution limit of the simulation and therefore many of 
the sink particles have an initial mass equal to the reso- 
lution limit of the simulation and all of the particles have 
an initial mass which is less than the initial Jeans mass. 
Therefore, the average mass in these isothermal simula- 
tions is unrelated to the initial Jeans mass (as discussed 
in MES06). For the simulations with dust heating, the 
initial sink masses at the beginning of the simulation 
are similar to those of the isothermal simulation, as Fig- 
ures [27] and [28] show (i.e., objects with high final masses 
form with low initial masses and early). However, as 
the simulation progresses and the temperature heats up, 
the initial sink masses become larger. It is important to 
note that even in this case, the initial sink masses are 
much smaller than the initial Jeans mass of our simula- 
tion, even though the average temperature has increased 
substantially from the initial 5K. This seems to indicate 
that in our simulations the initial Jeans mass does not 
predict the scale of fragmentation. 

Many of the sinks that form in the isothermal simu- 
lations are prevented from forming when dust heating 
is included. The gas that was destined to form these 
sinks instead gets eventually accreted by existing sinks, 
enabling them to reach large masses. As Table [4] shows, 
the final mass that ends up in sinks is comparable in the 
isothermal and dust simulations, but since the number of 
sinks is widely different, the masses are different as well. 

In Figures l2"2T[2"5l we also plot two observationally- 
derived mass functions. We use the analytic mass func- 



FlG. 27. — Initial versus final mass for Dust Heating simulation 
with N Rcn = 2. 
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Fig. 28. — Final sink mass as a function of order of formation 
(Sink Number) for Dust Heating simulation with iV gcn = 2. 

tions from ISalpeteri (fl955h and iChabrierl (f2003h (Disk 
IMF from Table 1) and normalize them to the maximum 
sink particle mass in each simulation at the final time. 
For the isothermal simulations, the mass function is very 
different from the two analytic mass functions. Too many 
low mass objects are formed given the mass of the most 
massive object in the simulation. However, the simula- 
tions with dust heating are able to form more massive 
stars. For the simulation with iVg Cn = 1, the mass func- 
tion appears to under-sample the analytic mass functions 
at low (M < IMq) and moderate masses (~ 1M Q ). The 
larger simulation (N gcn = 2) shows more promising re- 
sults. The mass function is better sampled. The slope 
at the high masses is very similar to the Salpeter slope. 
However under-sampling at low and intermediate masses 
is still present. 

At the same time, t = 2.5%, both simulations {N gcn — 
1 and 2), have a similar mass distribution with the same 
range of masses. As the simulation with N gea = 1 ad- 
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vances in time, its most massive objects gain mass and 
more low mass objects are formed. We cannot probe 
past t — 2.5iff in our iVg en = 2 simulation because the 
maximum allowable mass would be exceeded (see ij3.3[) . 
However, as we mention in §2.1[ we assume sink parti- 
cles do not undergo further fragmentation once they have 
formed. This may not be the case. If further fragmenta- 
tion does occur within some of our sink particles, these 
events may populate the low mass end of the distribution 
in Figure [231 

5. DISCUSSION 

The investigation performed in this work was a study 
of the effect of dust heating on the star formation process. 
This work is not intended to explain the complete star 
formation process. As we discuss in !|4j an obvious miss- 
ing factor in our work is t urbulence. Another impo rtant 
factor is magnetic fields (jPrice fe Batell2008l I2009D . In 
the following two sections we compare our work to other 
similar work and then discuss future improvements. 

5.1. Comparison to Previous Work 
In this section, we compare our work and results 



to those of 


Krumholz et al.l 


lOffher et al. 


(120091). The r 



tween our work and theirs is size. Our largest simula- 
tion (A'gon = 2) models a box with M = 671M Q and 
L = 0.984pc. It exceed s in size the largest simulation of 
IKrumholz et all (120071) which models 200M Q in a sphere 
of radius O.lpc or Offn er et al.l (j2009) with a box of mass 
185M Q and L = 0.65pc. Because of the differences 
in scale of our simulations, it is somewhat difficult to 
compare global properties since these other simulations 
do not form the large number of objects that we do. 
However, we all find that including the effect of radia- 
tive transfer drastically decreases the number of objects 
formed. (This can be understood from a Jeans mass 
argument. If the gas is hotter, in this case, due to ra- 
diation, then the amount of fragmentation will decrease, 
thereby reducing the number of objects formed.) So me of 
the oth er res ults that we fi nd ar e only hinted at in iBatd 
(|2009bf) and lOffner et al.l (|2009D . namely the change in 
the mass function and the correlation of accretion rate 
and mass. 

The highest num ber of objects formed in 
Krumholz et all (|2007l ) is 7. Therefore, they do 
not attem pt to produc e a m ass function. The same 
is true of Of fner et alJ (|2009T ) in which they form 15 
objects. IBatd (|2009bf ) forms 17 objects and produces 
a mass function which he compares to observed mass 
functions. However, with such a small number of objects 
and a maximum stellar mass less than 2M Q , it is difficult 
to draw conclusions. In our largest simulation with 
dust heating, we form 74 objects. Our mass function 
samples masses up to 20M©, and while accretion is still 
occurring in our simulation, we can see more conclusively 
that including dust heating encourages the formation 
of massive stars, while inhibiting the fragmentation 
that leads to an overabundance of low-mass objects in 
si mulations that do n ot include dust heating. 

lOffner et all (|2009l ) finds a slight trend of increasing 
accretion rate with mass. In IBatd (l2009al) (w hich uses 
similar simulation parameters as lBat c (2009b)), there is 
also a hint of a trend of accretion rate with mass. The 



trend in our simulation begins at ~ 2Mq and if objects 
larger than 10M Q are ignored then the trend is difficult 
to see in our smaller dust-heating simulation (see Fig. 
[T8|) which only formed 20 objects, sim ilar to the num- 
ber of objects formed in IBatd (|2009bl ) and lOffner et al.l 
(2009^. However, if we look at our larger simulation with 

74 objects, then the tren d is obvious (see Fig. \T§§. 

As d iscuss ed in CD IKrumholz et all (|2007h . IBatd 



(|2009bf ). and lOJher_eiaL| 420 09D all use FLD to calcu- 
late the dust temperature in their simula tions. Using this 
metho d in the optically thick regions of IKrumholz et al 



(120071) is pr obably valid. However, for the cases of iBate 
(|2009bh and lOffner et all (120091) . they study lower den- 
sity regions than IKrumholz et al.l (|2007f ) and it is unclear 
if the FLD approximation is still valid. All three also ig- 
nore the wavelength dependence of the dust opacity when 
calculating the radiation field. Our method also makes 
approximations, but of a different nature. We assume 
that the material is spherically distrib uted around the 
sink p articles. Based on th e figures in | Krumholz et alJ 
(pOOTT) . iBatel (l2009bf ). and lOffner et al.l (|2009l) this Is 
clearly not always the case. Both methods use approxi- 
mations, therefore, it is difficult to say which is a "better" 
method. Since there is currently no realistic method of 
using three-dimensional, wavelength-dependent radiative 
transfer, it may be the case that the best alternative is a 
combination of our two methods. At early stages when 
the gas is less dense and the density distribution around 
the sink particles is roughly spherical, our method may 
be more appropriate. However, as the density increases 
and disks begin to form around stars, the FLD method 
may be more appropriate. It is important to note that, 
despite the different radiative transfer methods used with 
a variety of approximatio ns, the main con c lusion of our 
work a nd th e works of | Krumholz et al.l (|2007f ). IBatd 
(|2009bh . and lOffner et alJ 1)20091 ) is the same: heating 
severely inhibits the fragmentation of the gas and pro- 
motes the formation of massive stars. 

One main difference between t he work oflBatd (12009b ) 
and our work (a n d the work of IKrumholz et all (2007) 



and lOffner et al.l ((2009)) is their a pproxi mation regard- 
ing the source of radiation. IBatd (|2009b) ignores much 
of the accretion luminosity and does not consider nuclear 
burning. The first assumption may not be valid during 
the early stages of star form ation because accretion may 
be at its highest le vel then. IKrumholz et alJ (|2007f) and 
lOffner et all (|2009l ) do not ignore the effect of protostel- 
lar heating, therefore their work do es not suffer f rom the 
problem of missing radiation as in Hate] (|2009bft . They 
both include a stellar model to calculate the accretion 
and intrinsic luminosity from the protostars in their sim- 
ulations. As discussed in the following section, ^5.21 we 
plan to use a similar approach using a stellar evolution 
model in our future work in order to improve our calcu- 
lation o f the stellar luminosity. 

Since IKrumholz etTaA (|2007h and lOffner et all (j2009l ) 
include a stellar model for their sink particles, we can 
compare the stellar properties from our simulations. If 
we compare our smaller simulation (N gcn — 1) to their 
simulations, we find that o ur luminosities are com para- 
ble. The ac c retion rate in IKrumholz et al.l (|2007l ) and 
lOffner et alJ ((2009) is found to be highly variable, simi- 
lar to our results. We can only qualitatively investigate 
the density distribution around the objects formed in the 
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simulations o flKrumholz et all (|2007t ). lBatd (|2009bf) . and 
lOffner et all (2009) via their simulation figures. From 
these it is not clear how to compare them with the 
spherically-averaged density distributions that we have 
calculated around our sink particles. However, they in- 
vestigate disk properties around some of their objects 
which we have not done. 

Another m a in difference between ou r wor k and that of 
IBatd (I2009H ). iKrumholz et alj (|2007t) . and lOffner et all 
(2009) is that our work does not include compressional 
heating or viscous heat i ng th at is included in the FLD 
method. lOffner et alj (|2009h find that compressional 
heating does not dominate the heating, but rather stellar 
heating is most important throughout most of their sim- 
ulation. However, before stars have formed, heating via 
viscous dissipation is dominant. Therefore, the effects 
that we have ignored are most likely to change how frag- 
mentation proceeds in the beginning of the simulation 
and then how the fir st sta r s form . Taken together, our 
work and the work of lBatd (|2009bf ) show the importance 
of different heating mechanisms throughout the star for- 
mation process. 

One minor differ ence between our work and that of 
lOffner et alj ([2009) is their inclusion of turbulent driv- 
ing throughout their simulation. Since our simulations 
are on different scales, it is not clear how this affects 
their results compared to ours. Another minor difference 
between our work and ot hers' is that we do no t allow 
sink pa rticle s to merge. [Krumholz et al.l (|2007ft . IBatd 
(|2009bf) . an d lOffner et all (|2009fl all allow this to occur; 
however for lBatd (|2009bf ) this never happens. Therefore, 
it is not clear if it is necessary to include this process in 
our simulations. 

5.2. Future Improvements 

Our method is an approximation of what we believe 
to be one of the most important effects in the very early 
stages of star formation, namely dust heating via young 
stars. We have attempted to model this stage as accu- 
rately as possible, yet there are areas which we believe 
can be improved in future work. 

(1) We do not account for all possible sources of radia- 
tion in our simulation. We have attempted to account for 
heating of the gas during the co llapse using the models of 
iWuchterl fc Tscharnuterl (|2003f ). However, there may be 
some extra contractional heating that we do not include 
when our sink particles have a mass less than O.OIMq. 

Besides including these effects, we can also improve 
the radiative transfer method discussed in this paper. 
Our method of calculating the dust temperature assumes 
spherical symmetry, yet images of young star forming re- 
gions show three-dimensional morphology. A more ad- 
vanced method of calculating the dust temperature, at 
the current level of computational power, is impractical. 
Even though we have already made approximations to 
reality, we have had to extrapolate in our dust temper- 
ature look-up table. Our extrapolation was necessitated 
by the fact that we could not know a priori the range 
of envelope profiles that would be created in our sim- 
ulation. This is discussed in more detail in W2A[ The 
extrapolation occurred in regions which we believe were 
fairly well-understood. We can address this in future 
work with an expanded interpolation table. 

(2) We currently assume that the dust and gas tem- 



perature are equal. Whil e this may be th e case for very 
dense regions studied by IKrumholz et all (|2007f) . this is 
not always the case for our model. We will address this 
issue in a forthcoming paper. 

(3) To calculate the luminosity, we have had to in- 
terpo late in a small table of mass and mass accretion 
rate (jWuchterl fc Tscharnuterl 12003). which introduced 
uncertainties in our luminosity calculation. Our method 
of calculating luminosity could be improved because we 
currently assume that luminosity varies smoothly with 
mass and mass accretion rate independent of past his- 
tory. We plan to address this in future work using an 
advanced stellar evolution code. 

(4) As mentioned before, fragmentation may be occur- 
ring within sink particles. We do not believe that this 
will strongly influence the temperature and thus the frag- 
mentation of the gas that has not yet formed into sinks. 
Consider the case of a high-mass sink that fragments 
(within the sink radius) into two equal or two highly 
dissimilar mass bodies. In the first case, the intrinsic lu- 
minosities from the two equal-mass stars would sum to 
less than the intrinsic luminosity of a single object with 
the same total mass because of the steep dependence of 
intrinsic luminosity on mass. However, the accretion lu- 
minosity would remain the same in both cases. As seen in 
Figure \13\ the accretion luminosity is a substantial con- 
tributor to the luminosity of objects with M > 10M©. 
Therefore, we do not expect the luminosity to change 
substantially in the first case. In the second case, the 
luminosity missing from the most massive sink due to 
fragmentation would be negligible. However, fragmenta- 
tion within a sink could affect the mass function of our 
simulation. This fragmentation could either decrease the 
mass of the most massive object in our simulation or it 
could populate the low-mass end of the mass function. 
This issue could be explored in future work by increas- 
ing the resolution of our simulation. 

6. SUMMARY AND CONCLUSION 

We have investigated the effect of the heating of dust 
via luminosity sources in a clustered star formation simu- 
lation. We compare the results of isothermal simulations 
to simulations that include dust heating. Including the 
effect of dust heating drastically reduces the number of 
objects formed (by more than an order of magnitude). 
We find that the density profiles of the envelopes sur- 
rounding the sinks/cores formed in our simulations with 
heating are comparable to those found around isolated, 
low-mass star-forming cores. This brings up the question 
of how similar density profiles can be formed in such dif- 
ferent accretion environments. Another interesting result 
of our simulations is that the accretion of mass onto the 
sinks /cores is found to be highly variable, in contrast 
to what is theorized for isolated, low-mass star-forming 
cores. We also find a strong correlation between the av- 
erage accretion rate and the final mass for objects with 
M > 2Mq. This fact may provide a clue to how massive 
stars form. We also analyze the final mass function of 
our simulations. We find that we are able to reproduce 
the results of MES06 for our isothermal simulations, i.e., 
a log-normal distribution centered at very low masses 
(~ O.IA/q) with no objects with masses greater than 
IMq. The mass functions produced by our simulations 
that include dust heating show that we are able to pro- 
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duce massive stars (M > 10M Q ). However, we do see a 
dearth of objects at low and intermediate masses. This 
may be due to the extreme heating by the dust and the 
lack of cooling physics. In our next paper we plan to 
relax the assumption of dust-gas collisional coupling at 
all densities. We will include the c omple te energetics al- 
gorithm described in lUrban et alj (|2009( ). which includes 
molecular cooling and cosmic-ray heating, in future sim- 
ulations similar to the ones discussed in this paper. 
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